function[x,output]=agent1solverho(x0,index,ewsb,waut);
global P beta n m S state agent gammagrid  Y scalefactor lowerbound rho table  ewdev statedev12 statedev13 statedev23
output=[];
objval = [];
options =optimset('LargeScale','off','Display','iter','MaxFunEvals',1000,'TolFun',0.00001);
[x,fval,exitflag]=fsolve(@objfun,x0,options);
output(6)=exitflag;
function f = objfun(x);
i=index(1);
j=index(2);
consu1 = x(1);
consu2 = x(2);

gammarfun=[x(3) x(4)];

wint=interpolation(scalefactor*ewsb, index,gammarfun);

output(1)=utility(rho,consu1)+beta*wint(1)/scalefactor;
output(2)=utility(rho,consu2)+beta*wint(2)/scalefactor;
output(3)=utility(rho,Y(j)-consu2-consu1)+beta*wint(3)/scalefactor;
f(1)=gammarfun(1) - consu2/consu1;
f(2)=gammarfun(2) - (Y(j)-consu2-consu1)/consu1;
f(3) = scalefactor*(utility(rho,S(j,1))+beta*P(j,:)*waut(:,1) - output(1));
f(4) = gammagrid(i,6)/gammagrid(i,5)-(Y(j)-consu2-consu1)/consu2;
end
end
